#include <iostream>
#include <vector>
#include <ctime>
#include <numeric>

#include "tensor.h"
#include "covariance.h"
#include "mtrand.h"
#include "dataset.h"
#include "pdf.h"
#include <vector>
#include <sstream>

#include "sdataset.h"
#include "spdf.h"

#include "TFile.h"
#include "TChain.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile2D.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TLorentzVector.h"
#include "TGraphErrors.h"
#include "TGraph2DErrors.h"
#include "make_list.h"

using namespace std;

//initialize random number generator
MTRand_int32 irand(time(0));

int main(int argc, char **argv)
{

  if(argc < 2+1)
    {
      cout<<"Usage: ./example_qcd input.root output.root"<<endl;
      throw;
    }
  
  
  //open rootfile
  TChain chain("Ntuple");
  chain.Add(argv[1]);

  //open output file
  TFile outfile(argv[2], "RECREATE");

  //declare variables
  const int maxsize=10;
  Int_t           size;
  Double_t         w;   //[size]
  Double_t         pt[maxsize];   //[size]
  Double_t         px[maxsize];   //[size]
  Double_t         py[maxsize];   //[size]
  Double_t         pz[maxsize];   //[size]
  Double_t         e[maxsize];   //[size]
  Double_t         m[maxsize];   //[size]
  Double_t         t1[maxsize];
  Double_t         t2[maxsize];
  Double_t         t3[maxsize];
  Double_t         t4[maxsize];


  //read the input file
  chain.SetBranchAddress("n", &size);
  chain.SetBranchAddress("weight", &w);
  chain.SetBranchAddress("pt", pt);
  chain.SetBranchAddress("px", px);
  chain.SetBranchAddress("py", py);
  chain.SetBranchAddress("pz", pz);
  chain.SetBranchAddress("e", e);
  chain.SetBranchAddress("m", m);
  chain.SetBranchAddress("tau1_1", t1);
  chain.SetBranchAddress("tau2_1", t2);
  chain.SetBranchAddress("tau3_1", t3);
  chain.SetBranchAddress("tau4_1", t4);

  //create some dataset
  //tau 21 of jet1 jet2 and deltaR
  const double PI=3.14159265359;

  DataSet DataSet_tau21_2jet(make_list<double>(0, 0, 1.0),
			     make_list<double>(1.0,1.0,2*PI),
			     make_list<double>(0.01,0.01,0.02));

  //DataSet for debugging
  DataSet DataSet_tau21_dR(make_list<double>(0, 1.0),
			   make_list<double>(1.0, 2*PI),
			   make_list<double>(0.01, 0.02));


  DataSet DataSet_tau21(make_list<double>(0),
			make_list<double>(1.0),
			make_list<double>(0.01));
  
  DataSet DataSet_tau21_jetmass(make_list<double>(0, 50.),
				make_list<double>(1.0, 500.),
				make_list<double>(0.01, 2.));

  DataSet DataSet_tau21_pt(make_list<double>(0, 50.),
			   make_list<double>(1.0, 500.),
			   make_list<double>(0.01, 2.));
  
  DataSet DataSet_jetmass_2jet(make_list<double>(50., 50., 1.0),
			       make_list<double>(500., 500., 2*PI),
			       make_list<double>(2, 2, 0.02));

  //declare some histograms
  map<string, TH1*> hist;
  hist["t21"]=new TH1D("t21", "tau 21", 20, 0, 1.0);

  hist["t21_dR"]=new TProfile2D("t21_dR", "tau 21 dR", 20, 0, 1.0,
				100, 1.0, 2*PI);

  hist["t21_MC"]=new TH1D("t21_MC", "tau 21", 20, 0, 1.0);
  hist["t21_MCerr"]=new TH1D("t21_MCerr", "tau 21", 20, 0, 1.0);
  hist["T21"]=new TH1D("T21", "T21 for two leading jet", 20, 0, 1.0);
  hist["T21_MC_jetmass"]=new TH1D("T21_jetmass", 
				  "T21 from jetmass PDF", 20, 0, 1.0);
  hist["T21_MCerr_jetmass"]=new TH1D("T21err_jetmass", 
				  "T21 err from jetmass PDF", 20, 0, 1.0);
  hist["T21_MC_pt"]=new TH1D("T21_pt", 
			     "T21 from pt PDF", 20, 0, 1.0);

  hist["T21_MC_corr_pt"]=new TH1D("T21_corr_pt", 
				  "T21 from pt PDF (corrected)", 20, 0, 1.0);
  hist["T21_MCerr_corr_pt"]=new TH1D("T21_corr_err_pt", 
				     "T21err from pt PDF (corrected)", 20, 0, 1.0);

  hist["T21_MCerr_pt"]=new TH1D("T21err_pt", 
			     "T21err from pt PDF", 20, 0, 1.0);

  hist["m1"]=new TH1D("m1", "mass for leading jet", 200, 50, 500);
  hist["m2"]=new TH1D("m2", "mass for second jet", 200, 50, 500);
  
  hist["T21_uncorrelated"]=new TH1D("T21_uncorrelated", 
				    "T21 for two leading jet", 20, 0, 1.0);
  

  int nmax=200000;

  //loop over events
  double total_weight=0;
  for(Long64_t i=0; i<chain.GetEntries(); i++)
    {
      
      if(i%10000==0)
      cout<<"\r\033[K"<<(i*10000/chain.GetEntries())*0.01<<"% trained..."
	  <<flush;


      chain.GetEntry(i);
      if(size<2) continue;
      
      //some sanity cuts
      if(t1[0] < 1e-5) continue;
      if(t1[1] < 1e-5) continue;
      if(t2[0] < 1e-5) continue;
      if(t2[1] < 1e-5) continue;
      if(pt[0] < 50.) continue;
      if(pt[1] < 50.) continue;
      if(m[0] < 50. || m[0] > 500.) continue;
      if(m[1] < 50. || m[1] > 500.) continue;

      TLorentzVector j1(px[0], py[0], pz[0], e[0]);
      TLorentzVector j2(px[1], py[1], pz[1], e[1]);

      double dR=j1.DeltaR(j2);
      if(dR>2*PI) continue;


      
      DataSet_tau21_jetmass.Fill(make_list<double>(t2[0]/t1[0], m[0]));
      DataSet_tau21_jetmass.Fill(make_list<double>(t2[1]/t1[1], m[1]));

      DataSet_tau21.Fill(make_list<double>(t2[0]/t1[0]));
      DataSet_tau21.Fill(make_list<double>(t2[1]/t1[1]));

      DataSet_tau21_pt.Fill(make_list<double>(t2[0]/t1[0], pt[0]));
      DataSet_tau21_pt.Fill(make_list<double>(t2[1]/t1[1], pt[1]));
      
      DataSet_tau21_2jet.Fill(make_list<double>(t2[0]/t1[0],
						t2[1]/t1[1],
						dR) );

      DataSet_tau21_dR.Fill(make_list<double>(t2[1]/t1[1],
					      dR) );

      DataSet_jetmass_2jet.Fill(make_list<double>(m[0], m[1], dR));

      total_weight+=w;
    }

  cout<<endl;

  //here compute come PDF
  PDF PDF_tau21_2jet=DataSet_tau21_2jet.ComputePDF(true);

  PDF PDF_tau21_dR=DataSet_tau21_dR.ComputePDF(true);
  PDF PDF_tau21_jetmass=DataSet_tau21_jetmass.ComputePDF(true);
  PDF PDF_tau21_pt=DataSet_tau21_pt.ComputePDF(true);
  PDF PDF_tau21=DataSet_tau21.ComputePDF(true);
  PDF PDF_jetmass_2jet=DataSet_jetmass_2jet.ComputePDF(true);

  //get MC_Generator for computation
  PDF::MC_Generator MC_tau21_2jet= PDF_tau21_2jet.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_jetmass_2jet= PDF_jetmass_2jet.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  //initialize some generator
  PDF::MC_Generator MC_tau21_jetmass
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::INPUT), 
     &PDF_tau21_jetmass);

  PDF::MC_Generator MC_tau21_pt
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::INPUT),
     &PDF_tau21_pt);
  
  PDF::MC_Generator MC_tau21
    (make_list<PDF::MC_Flag>(PDF::OUTPUT),
     &PDF_tau21);

  PDF::MC_Generator MC_tau21_dR
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::OUTPUT),
     &PDF_tau21_dR);


  double weight_cut=total_weight*0.01/chain.GetEntries();
  //loop over events a second time 
  for(Long64_t i=0; i<chain.GetEntries(); i++)
    {
      if(i%10000==0)
      cout<<"\r\033[K"<<(i*10000/chain.GetEntries())*0.01<<"% computed..."
	  <<flush;

      chain.GetEntry(i);

      if(w<weight_cut) continue;

      //////////////////
      //// **********
      //// Set Weights to 1
      /////////////////
      
      if(size<2) continue;
      
      //some sanity cuts
      if(t1[0] < 1e-5) continue;
      if(t1[1] < 1e-5) continue;
      if(t2[0] < 1e-5) continue;
      if(t2[1] < 1e-5) continue;
      if(pt[0] < 50.) continue;
      if(pt[1] < 50.) continue;
      if(m[0] < 50. || m[0] > 500.) continue;
      if(m[1] < 50. || m[1] > 500.) continue;

      TLorentzVector j1(px[0], py[0], pz[0], e[0]);
      TLorentzVector j2(px[1], py[1], pz[1], e[1]);

      double dR=j1.DeltaR(j2);
      if(dR>2*PI) continue;

      double t21_1=t2[0]/t1[0];
      double t21_2=t2[1]/t1[1];

      hist["t21"]->Fill(t21_1);
      hist["t21"]->Fill(t21_2);

      hist["T21"]->Fill(sqrt(t21_1*t21_2), w);

      int nGenerate=200;

      //generate some toy MC
      vector<PDF::MC_Event> MC_jet1=
	MC_tau21_jetmass.GenerateFast(make_list<double>(m[0]), nGenerate);

      vector<PDF::MC_Event> MC_jet2=
	MC_tau21_jetmass.GenerateFast(make_list<double>(m[1]), nGenerate);

      //since we are doing something complicated
      //we need to rescale all weight to 1
      double rescale=0;
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	rescale+= MC_jet1[j].weight* MC_jet2[j].weight;

      rescale=1.0/rescale;      
      
      //now loop over the MC_event
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  double temp_t21_1=MC_jet1[j].value[0];
	  double temp_t21_2=MC_jet2[j].value[0];
	  double temp_T21=sqrt(temp_t21_1*temp_t21_2);

	  double temp_weight=MC_jet1[j].weight*MC_jet2[j].weight*rescale;
	  if(temp_weight==0) continue;
	  double temp_err1=MC_jet1[j].err/MC_jet1[j].weight;
	  double temp_err2=MC_jet2[j].err/MC_jet2[j].weight;
	  
	  double temp_err=temp_weight*sqrt(temp_err1*temp_err1+
					   temp_err2*temp_err2);

	  hist["T21_MC_jetmass"]->Fill(temp_T21, w*temp_weight);
	  hist["T21_MCerr_jetmass"]->Fill
	    (temp_T21, w*temp_err);
	}

      //do the same thing for pt
      MC_jet1=
	MC_tau21_pt.GenerateGrad(make_list<double>(pt[0]), nGenerate);

      MC_jet2=
	MC_tau21_pt.GenerateGrad(make_list<double>(pt[1]), nGenerate);


      //get the correlation
      double correlation=0;

      vector<PDF::MC_Event> MC_corr=
	MC_tau21_2jet.GenerateFull(make_list<double>(j1.DeltaR(j2)));
      
      double mean_xy=0;
      double mean_x=0;
      double mean_y=0;
      double sum_weight=0;

      for(unsigned int j=0; j<MC_corr.size(); ++j)
	{
	  double temp_x=MC_corr[j].value[0];
	  double temp_y=MC_corr[j].value[1];
	  double temp_weight=MC_corr[j].weight;
	  sum_weight+=temp_weight;
	  mean_xy+=temp_x*temp_y*temp_weight;
	  mean_x+=temp_x*temp_weight;
	  mean_y+=temp_y*temp_weight;
	}

      sum_weight=1.0/sum_weight;
      correlation=(mean_xy*sum_weight - mean_x*mean_y*sum_weight*sum_weight);
      

      //since we are doing something complicated
      //we need to rescale all weight to 1
      rescale=0;
      double rescale_reweight=0;
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  rescale+= MC_jet1[j].weight* MC_jet2[j].weight;
	  //get the correction
	  rescale_reweight +=
	    MC_jet1[j].weight* MC_jet2[j].weight*
	    (1+ correlation*
	     MC_jet1[j].GetExtra()[0] * MC_jet2[j].GetExtra()[0]);
	}


      rescale=1.0/rescale;      
      rescale_reweight=1.0/rescale_reweight;
      
      //now loop over the MC_event
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  double temp_t21_1=MC_jet1[j].value[0];
	  double temp_t21_2=MC_jet2[j].value[0];
	  double temp_T21=sqrt(temp_t21_1*temp_t21_2);

	  double temp_weight=MC_jet1[j].weight*MC_jet2[j].weight*rescale;
	  if(temp_weight==0) continue;
	  double temp_err1=MC_jet1[j].err/MC_jet1[j].weight;
	  double temp_err2=MC_jet2[j].err/MC_jet2[j].weight;
	  
	  double temp_err=temp_weight*sqrt(temp_err1*temp_err1+
					   temp_err2*temp_err2);

	  //get the correction
	  double reweight=(1+ correlation*
			   MC_jet1[j].GetExtra()[0] * 
			   MC_jet2[j].GetExtra()[0]);

	  double temp_weight_reweight=
	    reweight*MC_jet1[j].weight*MC_jet2[j].weight*rescale_reweight;

	  double temp_err_reweight=temp_weight_reweight*
	    sqrt(temp_err1*temp_err1+
		 temp_err2*temp_err2);

	  hist["T21_MC_pt"]->Fill(temp_T21, w*temp_weight);
	  hist["T21_MC_corr_pt"]->Fill(temp_T21, w*temp_weight_reweight);

	  hist["T21_MCerr_corr_pt"]->Fill(temp_T21, w*temp_err_reweight);

	  hist["T21_MCerr_pt"]->Fill
	    (temp_T21, w*temp_err);
	}

      hist["m1"]->Fill(m[0], w);
      hist["m2"]->Fill(m[1], w);
    }

  cout<<endl;


  //generate t21 MC
  vector<PDF::MC_Event> MC_evt=
    MC_tau21.GenerateFast(100000);  
  

  //loop over MC Events
  for(unsigned int i=0; i<MC_evt.size(); ++i)
    {
      hist["t21_MC"]->Fill(MC_evt[i].value[0], MC_evt[i].weight);
      hist["t21_MCerr"]->Fill(MC_evt[i].value[0], MC_evt[i].err);
    }

  //generate t21_dR MC
  vector<PDF::MC_Event> MC_dR_evt=
    MC_tau21_dR.GenerateFull(make_list<double>());
  
  for(unsigned int i=0; i<MC_dR_evt.size(); ++i)
    {
      ((TProfile2D*) hist["t21_dR"] )->
	Fill(MC_dR_evt[i].value[0], MC_dR_evt[i].value[1], 
	     MC_dR_evt[i].weight);
    }
  


  //generate MC for correlation studies
  //first grab some more histograms
  hist["t21_corr"]=new TProfile("t21_corr", "#tau_{21} Correlation", 100, 1.0, 2*PI);
  hist["t21_info"]=new TProfile("t21_info", "#tau_{21} Total Information", 100, 1.0, 2*PI);


  hist["jetmass_corr"]=new TProfile("jetmass_corr", "Jet mass Correlation", 100, 1.0, 2*PI);
  hist["jetmass_info"]=new TProfile("jetmass_info", "Jet mass Total Information", 100, 1.0, 2*PI);


  for(float i=1.2; i<2*PI; i+= ((2*PI-1.2)/100.))
    {
      vector<PDF::MC_Event> MC_t21_evt=
	MC_tau21_2jet.GenerateFull(make_list<double>(i));
      
      //DEBUG
      cout<<"MC size "<<MC_t21_evt.size()<<endl;

      //this is really bad code for now
      //heavily depends on binning schemes used before
      //do not use for other things
      //loop over MC events
      //cout the size and make sure it is a square root
      int sqrt_size=sqrt(MC_t21_evt.size());
      int sqrt_size_sq=sqrt_size*sqrt_size;
      
      assert(MC_t21_evt.size()==sqrt_size_sq);

      //now get projected PDF
      valarray<double> PDFx(sqrt_size);
      valarray<double> PDFy(sqrt_size);
      
      for(unsigned int j=0; j<sqrt_size; ++j)
      for(unsigned int k=0; k<sqrt_size; ++k)
	{
	  PDFx[j]+=MC_t21_evt[j+k*sqrt_size].weight;
	  PDFy[k]+=MC_t21_evt[j+k*sqrt_size].weight;	  
	}
      
      double sumx=PDFx.sum();
      double sumy=PDFy.sum();

      //get the new bin volume for normalization
      double new_bin_vol=DataSet_tau21_2jet.width[1];

      cout<<"volume is : "<<new_bin_vol<<endl;

      //rescale the PDF
      if(sumx>0)
	{
	  sumx=1.0/(sumx*new_bin_vol);
	  PDFx*=sumx;
	}
      if(sumy>0)
	{
	  sumy=1.0/(sumy*new_bin_vol);
	  PDFy*=sumy;
	}

      
      //first we compute the correlation
      //correlation = (mean xy - mean x*meany)/ (sigmax* sigmay)
      double mean_xy=0;
      double mean_x=0;
      double mean_y=0;
      double mean_xx=0;
      double mean_yy=0;
      double sum_weight=0;
      
      for(unsigned int j=0; j<MC_t21_evt.size(); ++j)
	{
	  double temp_x=MC_t21_evt[j].value[0];
	  double temp_y=MC_t21_evt[j].value[1];
	  double temp_weight=MC_t21_evt[j].weight;
	  sum_weight+=temp_weight;
	  mean_xy+=temp_x*temp_y*temp_weight;
	  mean_x+=temp_x*temp_weight;
	  mean_y+=temp_y*temp_weight;
	  mean_xx+=temp_x*temp_x*temp_weight;
	  mean_yy+=temp_y*temp_y*temp_weight;
	}

      //make sure sum_weight > 0
      //to fill histogram
      if(sum_weight>0)
	{
	  //DEBUG
	  cout<<"sum_weight is : "<<sum_weight<<endl;
	  sum_weight=1.0/sum_weight;
	  mean_xy*=sum_weight;
	  mean_x*=sum_weight;
	  mean_y*=sum_weight;
	  mean_xx*=sum_weight;
	  mean_xy*=sum_weight;
	 
	  //DEBUG
	  cout<<"i,x,y,xx,yy,xy: "<<i<<","<<mean_x<<","<<mean_y<<","
	      <<mean_xx<<","<<mean_yy<<","<<mean_xy<<endl;

	  cout<<"xrms, yrms: "<<sqrt(mean_xx - mean_x*mean_x)<<","
	      <<sqrt(mean_y - mean_y*mean_y)<<endl;

	  cout<<"corr: "<<	     (mean_xy - mean_x*mean_y)/
	     sqrt(
		  (mean_xx - mean_x*mean_x)*
		  (mean_yy - mean_y*mean_y)
		  )
	      <<endl;
 
	  hist["t21_corr"]->Fill
	    (i, 
	     (mean_xy - mean_x*mean_y)/
	     sqrt(
		  (mean_xx - mean_x*mean_x)*
		  (mean_yy - mean_y*mean_y)
		  )
	     );
	  
	}


      //now compute total mutual information
      //integral rho(x,y) *log( rho(x,y)/(rho(x)*rho(y)) );
      double volume=0.01*0.01;
      double integral=0;
      
      for(unsigned int j=0; j<sqrt_size; ++j)
	for(unsigned int k=0; k<sqrt_size; ++k)
	  {
	    double rho_xy=MC_t21_evt[j+k*sqrt_size].weight;
	    double rho_x=PDFx[j];
	    double rho_y=PDFy[k];
	    
	    if(rho_xy >0 && rho_x >0 && rho_y >0)
	      integral-=
		volume*rho_xy*log(rho_xy/(rho_x*rho_y)); 
	  }
      
      //we have our total information
      hist["t21_info"]->Fill(i, integral);
      
      
      
      


    };



  //jetmass
  for(float i=1.2; i<2*PI; i+= ((2*PI-1.2)/100.))
    {
      vector<PDF::MC_Event> MC_t21_evt=
	MC_jetmass_2jet.GenerateFull(make_list<double>(i));
      
      //DEBUG
      cout<<"MC size "<<MC_t21_evt.size()<<endl;

      //this is really bad code for now
      //heavily depends on binning schemes used before
      //do not use for other things
      //loop over MC events
      //cout the size and make sure it is a square root
      int sqrt_size=sqrt(MC_t21_evt.size());
      int sqrt_size_sq=sqrt_size*sqrt_size;
      
      assert(MC_t21_evt.size()==sqrt_size_sq);

      //now get projected PDF
      valarray<double> PDFx(sqrt_size);
      valarray<double> PDFy(sqrt_size);
      
      for(unsigned int j=0; j<sqrt_size; ++j)
      for(unsigned int k=0; k<sqrt_size; ++k)
	{
	  PDFx[j]+=MC_t21_evt[j+k*sqrt_size].weight;
	  PDFy[k]+=MC_t21_evt[j+k*sqrt_size].weight;
	}
      
      double sumx=PDFx.sum();
      double sumy=PDFy.sum();

      //get the new bin volume for normalization
      double new_bin_vol=DataSet_jetmass_2jet.width[1];

      //rescale the PDF
      if(sumx>0)
	{
	  sumx=1.0/(sumx*new_bin_vol);
	  PDFx*=sumx;
	}
      if(sumy>0)
	{
	  sumy=1.0/(sumy*new_bin_vol);
	  PDFy*=sumy;
	}

      
      //first we compute the correlation
      //correlation = (mean xy - mean x*meany)/ (sigmax* sigmay)
      double mean_xy=0;
      double mean_x=0;
      double mean_y=0;
      double mean_xx=0;
      double mean_yy=0;
      double sum_weight=0;
      
      for(unsigned int j=0; j<MC_t21_evt.size(); ++j)
	{
	  double temp_x=MC_t21_evt[j].value[0];
	  double temp_y=MC_t21_evt[j].value[1];
	  double temp_weight=MC_t21_evt[j].weight;
	  sum_weight+=temp_weight;
	  mean_xy+=temp_x*temp_y*temp_weight;
	  mean_x+=temp_x*temp_weight;
	  mean_y+=temp_y*temp_weight;
	  mean_xx+=temp_x*temp_x*temp_weight;
	  mean_yy+=temp_y*temp_y*temp_weight;
	}

      //make sure sum_weight > 0
      //to fill histogram
      if(sum_weight>0)
	{
	  //DEBUG
	  cout<<"sum_weight is : "<<sum_weight<<endl;
	  sum_weight=1.0/sum_weight;
	  mean_xy*=sum_weight;
	  mean_x*=sum_weight;
	  mean_y*=sum_weight;
	  mean_xx*=sum_weight;
	  mean_xy*=sum_weight;
	 
	  //DEBUG
	  cout<<"i,x,y,xx,yy,xy: "<<i<<","<<mean_x<<","<<mean_y<<","
	      <<mean_xx<<","<<mean_yy<<","<<mean_xy<<endl;
 
	  hist["jetmass_corr"]->Fill
	    (i, 
	     (mean_xy - mean_x*mean_y)/
	     sqrt(
		  (mean_xx - mean_x*mean_x)*
		  (mean_yy - mean_y*mean_y)
		  )
	     );
	  
	}


      //now compute total mutual information
      //integral rho(x,y) *log( rho(x,y)/(rho(x)*rho(y)) );
      double volume=(500.-50.)*(500.-50.)/4.;
      double integral=0;
      
      for(unsigned int j=0; j<sqrt_size; ++j)
	for(unsigned int k=0; k<sqrt_size; ++k)
	  {
	    double rho_xy=MC_t21_evt[j+k*sqrt_size].weight;
	    double rho_x=PDFx[j];
	    double rho_y=PDFy[k];
	    
	    if(rho_xy >0 && rho_x >0 && rho_y >0)
	      integral-=
		volume*rho_xy*log(rho_xy/(rho_x*rho_y)); 
	  }
      
      //we have our total information
      hist["jetmass_info"]->Fill(i, integral);
      
      
      
      


    };



  //write all histograms
  map<string, TH1*>::iterator it;
  for(it=hist.begin(); it!=hist.end(); ++it)
    {
      //see if the name has MC but no MCerr
      if(it->first.find("MC")!=string::npos &&
	 it->first.find("MCerr")==string::npos)
	{
	  //try to see if the corresponding MCerr histogram exists
	  string str=it->first;
	  str.replace(str.find("MC"), 2, "MCerr");
	  TH1* ch=hist[str];
	  if(ch!=NULL)
	    {
	      //found a corresponding histogram
	      for(int i=1; i<=it->second->GetNbinsX(); ++i)
		{
		  it->second->SetBinError
		    (i, ch->GetBinContent(i));
		}

	    }
	}

      it->second->Write();
    }

  outfile.Close();
  cout<<"Program ends."<<endl;
  return 1;
}
  
